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The DAMA/Nal and DAMA/LIBRA annual modulation data, which may be in- 
terpreted as a signal for the existence of weakly interacting dark matter (WIMPs) 
in our galactic halo, are examined in light of null results from other experiments: 
CDMS, XENON10, CRESST I, CoGeNT, TEXONO, and Super-Kamiokande (Su- 
perK). We use the energy spectrum of the combined DAMA modulation data given 
in 36 bins, and include the effect of channeling. Several statistical tools are imple- 
mented in our study: likelihood ratio with a global fit and with raster scans in the 
WIMP mass and goodness-of-fit (g.o.f.). These approaches allow us to differentiate 
between the preferred (global best fit) and allowed (g.o.f.) parameter regions. It is 
hard to find WIMP masses and couplings consistent with all existing data sets; the 
surviving regions of parameter space are found here. For spin-independent (SI) inter- 
actions, the best fit DAMA regions are ruled out to the 3a C.L., even with channeling 
taken into account. However, for WIMP masses of ~8 GeV some parameters outside 
these regions still yield a moderately reasonable fit to the DAMA data and are com- 
patible with all 90% C.L. upper limits from negative searches, when channeling is 
included. For spin-dependent (SD) interactions with proton-only couplings, a range 
of masses below 10 GeV is compatible with DAMA and other experiments, with 
and without channeling, when SuperK indirect detection constraints are included; 
without the SuperK constraints, masses as high as ~20 GeV are compatible. For SD 
neutron-only couplings we find no parameters compatible with all the experiments. 
Mixed SD couplings are examined: e.g. ~8 GeV mass WIMPs with a n = ±a p are 
found to be consistent with all experiments. In short, there are surviving regions at 
low mass for both SI and SD interactions; if indirect detection limits are relaxed, 
some SD proton-only couplings at higher masses also survive. 
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I. INTRODUCTION 



Among the best motivated candidates for dark matter are Weakly Interacting Massive 
Particles (WIMPs). Direct searches for dark matter WIMPs aim at detecting the scattering 
of WIMPs off of nuclei in a low-background detector. These experiments measure the energy 
of the recoiling nucleus, and are sensitive to a signal above a detector-dependent energy 
threshold [lj]. 

The discovery of an annual modulation by the DAMA/Nal experiment 0], confirmed 
by the new experiment DAMA/LIBRA Q of the same collaboration, is the only positive 
signal seen in any dark matter search. The DAMA collaboration has found an annual 
modulation in its data compatible with the signal expected from dark matter particles bound 
to our Galactic Halo 0, H|. However, other direct detection experiments, e.g. CDMS (1,01, 
CoGeNT (slB^OUPP [H, 13, CRESST 0, [3, KIMS Q, TEXONO [H, 3, and 
XENON10 [171. Il8j| have not found any signal from WIMPs. It has been difficult to reconcile 
a WIMP signal in DAMA with the other negative results [19], for the case of canonical 
WIMP masses ~ lOOGeV with standard weak interactions motivated by Supersymmetry 
(SUSY). Yet, for light WIMPs, such compatibility was possible. Papers written prior to the 
latest DAMA/LIBRA results found regions in cross-section and WIMP mass that reconciled 
all null results with DAMA's positive si gnal , even assuming a standard halo model: WIMPs 
with spin independent interactions 



20, 21 



22 



in the mass range 5-9 GeV, and WIMPs with 
in the mass range 5-13 GeV were found to be compatible 
Ref. [H, 21 1 studied not only the case of the standard 



spin dependent interactions 
with all existing data. In addition 
halo model, but also models with dark matter streams; they found that nonstandard velocity 
distributions could reconcile all the existing data sets including DAMA. Other dark matter 
explanations invoke inelastic scattering of WIMPs (23| . 

Since these studies of a few years ago, DAMA/LIBRA has reported new experimental 
results. First, of particular interest is the new presentation of the modulation data in 36 
separate energy bins, which allows much more detailed investigations of the data. Second, 
in the past year another possible complication has come to light. The DAMA collaboration 
has pointed out that "channeling" may affect the interpretation of their data. In general 
only a fraction (known as the quenching factor) of the recoil energy deposited by a WIMP 
is transferred to electrons and is converted into useful signal in the DAMA detector (e.g. 
ionization or scintillation); the remainder is converted into phonons and heat and goes un- 
detected. Hence measured energies must be corrected for this behavior to obtain the proper 
recoil energy, by dividing by this quenching factor. The DAMA detector is composed of Nal, 
with quenching factors Qno. = 0.3 and Qi = 0.09. Yet, recently, the DAMA collaboration 
pointed out that some fraction of the nuclear recoils travel in straight paths along the crystal 
in such a way as to lose very little energy to other atoms and to heat, giving nearly all their 
energy to electrons; i.e., for these cases the measured energy corresponds very nearly to the 
energy deposited by the WIMP and thus have Q ~ 1. As a consequence, the detector is 
sensitive to lower mass WIMPs than previously thought. The DAMA threshold is 2 keVee 
(electron equivalent). Using the above quenching factors this was thought to correspond to 
7 keV recoil energy for sodium and 22 keV recoil energy for iodine; yet, due to channeling 
the actual threshold for some of the events can be as low as 2 keV. The new channeling 
studies revived the possiblity of DAMA's compatibility with other data sets, particularly at 
low masses. 

have been 



Light neutralinos as WIMPs with masses as low as 2 GeV 24] or 6 GeV [25l . |26 
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considered, even with the cross sections we find necessary for WIMPs to be compatible with 
the DAMA signal and all other negative searches (for spin independent interactions) [26]. 
In any event, in this paper we proceed in a purely phenomeno logical way in choosing the 
WIMP mass and cross section, with either spin-independent or spin-dependent cross sec- 
tions. We do not attempt to provide an elementary particle model to support the values 
of masses and cross sections. As justification of our approach, let us recall that there is no 
proven particle theory of dark matter. The candidates we are considering are stable neutral 
particles which have very small cross sections with nucleons. Regarding their production in 
accelerators, they would escape from the detectors without interacting. Unless there is a 
concrete specific model relating our neutral candidate to other charged particles (which can, 
if fact, be observed) there is no way such particles could be found in accelerators. The usual 
signature searched for in accelerators, for example at LEP, Tevatron or LHC, is the emission 
of a charged particle related to the neutral particle in question. For example, searching 
for "neutralinos" one puts bounds on one of its cousins, a "chargino" , or another relative, a 
"slepton" . Without a detailed model there are no accelerator bounds on neutral dark matter 
candidates. 

An alternative way to search for WIMP dark matter of relevance to this paper is via 
indirect detection of WIMP annihilation in the Sun. When WIMPs pass through a large 
celestial body, such as the Sun or the Earth 27|, |28( , interactions can lead to gravitational 
capture if enough energy is lost in the collision to fall below the escape velocity. Captured 
WIMPs soon fall to the body's core and eventually annihilate with other WIMPs. These 
annihilations lead to high-energy neutrinos that can be observed by Earth-based detectors 
such as Super-Kamiokande [29[ (SuperK), which currently provides the tightest indirect de- 
tection bounds at light WIMP masses, AMANDA IceCube Q and ANTARES [33 



The annihilation rate depends on the capture rate of WIMPs, which is in turn determined 
by the WIMP scattering cross section off nuclei in the celestial body. While the Earth is 
predominantly composed of spinless nuclei, the Sun is mostly made of hydrogen, which has 
spin. Thus the spin-dependent cross section of WIMPs off nucleons can be probed by mea- 
suring the annihilation signals from WIMP annihilation in the Sun. Other indirect detection 
methods search for WIMPs that annihilate in the Galactic Halo or near the Galactic Center 
where they produce neutrinos, positrons, or antiprotons that may be seen in detectors on the 



Earth [3J, [35(. Of course, these rely upon additional physics (WIMP annhilation) beyond 
what direct detection experiments rely on (WIMP scattering). Thermal WIMPs can only 
achieve the correct relic density today if their annihilation cross sections are fixed to be the 
value we use when computing the SuperK limits; yet one can imagine alternatives in which 
case the SuperK bounds should not be used. 

In this paper, we reconsider the issue of compatibility of the new DAMA/Nal and 
DAMA/LIBRA (hereafter DAMA) results with all other experimental results. We use the 
36 bins of modulation data and take into account the channeling effect. We restrict our 
studies to the standard isothermal model for the dark matter halo. We carefully investigate 
both spin-independent as well as spin-dependent couplings. Our procedure is the following. 
Assuming the standard dark halo model, we find the viable region by making sure that we 
produce the correct amplitudes for the DAMA modulation. We then compare this region 
with the most stringent current experimental constraints from negative searches. We use 
three different statistical methods to find the DAMA region: the likelihood method with (i) 
a global fit and (ii) raster scans in the WIMP mass, and (iii) goodness-of-fit. These methods 
and the different information derived from them are described in detail. 
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II. DARK MATTER DETECTION 



WIMP direct detection experiments seek to measure the energy deposited when a WIMP 
interacts with a nucleus in the detector |36j . If a WIMP of mass m scatters elastically 
from a nucleus of mass M, it will deposit a recoil energy E = (fi 2 v 2 /M)(l — cos 6), where 
fi = mM/ (m + M) is the reduced mass, v is the speed of the WIMP relative to the nucleus, 
and 9 is the scattering angle in the center of mass frame. The differential recoil rate per 
unit detector mass for a WIMP mass m, typically given in units of cpd/kg/keV (where cpd 
is counts per day), can be written as: 

dR a(q) . < 



where q = y2ME is the nucleus recoil momentum, a(q) is the WIMP-nucleus cross-section, 
p is the local WIMP density, and information about the WIMP velocity distribution is 
encoded into the mean inverse speed rj(E,t), 

v(E,t)= [ (2) 



u>v n 



Here 



represents the minimum WIMP velocity that can result in a recoil energy E and /(u, t) is 
the (time-dependent) distribution of WIMP velocities u relative to the detector. 

To determine the number of expected recoils for a given experiment and WIMP mass, we 
integrate Eqn. ([1]) over the nucleus recoil energy to find the recoil rate R per unit detector 
mass: 

rE 2 /Q 

R(t)= dEe(QE)-^a(q) V (E,t). (4) 

JE 1 /Q 

e(QE) is the (energy dependent) efficiency of the experiment, due, e.g., to data cuts designed 
to reduce backgrounds. Q is the quenching factor relating the observed energy E det (in 
some cases referred to as the electron-equivalent energy) with the actual recoil energy E rec \ 
Edet — QE rec (explained in more detail below). The energy range between E\ and E 2 is 
that of observed energies for some data bin of the detector (where experiments often bin 
observed recoils by energy). The quenching factor Q depends on the nuclear target and the 
characteristics of the detector. 

The recoiling nucleus will transfer its energy to either electrons, which may be observed 
as e.g. ionization or scintillation in the detector, or to other nuclei, producing phonons and 
heat. Experiments that do not measure the phonons/heat can only directly measure the 
fraction of energy Q that goes into the channel that is observed (such as scintillation); we 
refer to this as the observed energy. Observed energies for experiments that measure only 
from the electron channel will be quoted in electron-equivalent energies (keVee). Note some 
detectors can calibrate their energy scales to Edet = E rec , in which case the quenching factor 
can be ignored in the above formulations (such calibrations may directly or indirectly involve 
a quenching factor). 
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Real experimental apparatus cannot determine the event energies with perfect precision. 
If the expected amount of energy in the channel an experiment observes due to a nuclear 
recoil is E', the measured energy will be distributed about E'. The previous formulas apply 
only for a perfect energy resolution and must be corrected for this finite resolution. The 
observed rate over a measured energy range E\-E<i, Eqn. (jlj), should be replaced with 

R(t) = [°° dE e(QE) $(QE, E u E 2 ) a{q) V (E, t). (5) 

Jo 2m V 

Here, = QE, E±, E2) is a response function corresponding to the fraction of events 

with an expected observed energy E' = QE (i.e. the energy that would be observed with 
perfect energy resolution) that will actually be measured between E\ and E 2 ; recall E in 
the above equation is the recoil energy, of which only a fraction Q is converted into a form 
that is potentially observable. If the measured energies are normally distributed about E' 
with a (typically energy dependent) standard deviation cr(E'), 



$(E',E 1 ,E 2 ) = ± 



erf I — % 1 — erf 1 



V2a(E')J \V2ct(E' 



(6) 



In the limit cr(E') — > (perfect energy resolution), Eqn. (jSJ) reduces to Eqn. (TjJ. 
For detectors with multiple elements and/or isotopes, the total rate is given by: 



Rtat(t)=J2fi R i(t) (7) 

i 

where /j is the mass fraction and Ri is the rate (Eqn. for element /isotope i. 
The expected number of recoils observed by a detector is given by: 

Nrec = M det TR (8) 

where M^et is the detector mass and T is the exposure time. 



A. Cross-Section 

The a(q) cross-section term in Eqns. (pQ) &; @ is an effective cross-section for scatters 
with a momentum exchange q. The momentum exchange dependence appears in form factors 
that arise from the finite size of the nucleus. The total scattering cross-section generally has 
contributions from spin-independent (SI) and spin-dependent (SD) couplings, with 

a = cr SI + o- SD ; (9) 

these two cross-sections are described below. 

Spin-independent (SI). For spin-independent WIMP interactions, we make the usual as- 
sumption [l| that the cross section a scales with the square of the nucleus atomic number 
A and is given by 

a = a \E(E)\ 2 (10) 

where <7o is the zero-momentum WIMP-nuclear cross-section and F(E) is a nuclear form 
factor, normalized to F(0) = 1. For purely scalar interactions, 

4/y 2 

a 0iS1 = -^[Zf p +(A-Z)f n } 2 . (11) 
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Here Z is the number of protons, A — Z is the number of neutrons, and / p and f n are the 
WIMP couplings to the proton and nucleon, respectively. In most instances, f n ~ f p ; the 
WIMP-nucleus cross-section can then be given in terms of the WIMP-proton cross-section 
as a result of Eqn. (FiT|) : 

' JL 

where the /x p is the proton- WIMP reduced mass, and A is the atomic mass of the target 
nucleus. In this model, for a given WIMP mass, ct p si is the only free parameter. 
For the nuclear form factor we use the conventional Helmi form 37]. 



^o.si 



CT 



p,SI 



A 



(12) 



F(E) = 3e 



-q 2 s 2 /2 



sin(gr) — qr cos(qr) 
(qr) 3 



(13) 



with s = 0.9 fm, r is an effective nuclear radius as described in Ref. [37]] , and q = \p2ME. 
Spin- dependent (SD). The generic form for the spin-dependent WIMP-nucleus cross- 



section includes two couplings 
coupling a n , 



[38l | . the WIMP-proton coupling a p and the WIMP-neutron 



<?SD(q) 



2J+1 



[a^Sppiq) + a p a n S pn (q) + a 2 n S nn (q)\ 



(14) 



Here, the quantities Op and a, 
in units of 2y/2G F H |j 



are actually the axial four-fermion WIMP-nucleon couplings 



P n(q) 



411 ] . The nuclear structure functions S pp (q), S nn (q), and S, 
are functions of the exchange momentum q and are specific to each nucleus. We take 
the structure functions for Aluminum from Ref. 42l 
Ref. 43j; and for Silicon and Germanium from Ref. 
functions may be found in the review of Ref. d 



for Sodium, Iodine, and Xenon from 
I44J1. These and additional structure 



B. Velocity Distribution 

The dark matter halo is likely to have a smooth, well mixed (virialized) component that 
will lead to detectable recoils in direct detection experiments; the simplest model of this 
component is the Standard Halo Model (SHM) In addition to this smooth component, 
the galaxy contains structure from the galaxy formation process that has not yet virialized. 
This structure includes, e.g., tidal streams of dwarf galaxies in t he p rocess of being absorbed 
by the Milky Way, such as the Sagittarius dwarf galaxy j46|, |47], S, 0, 50]. Any such 



structure located about the Solar System will also produce scattering events. 

For each of these components, the smooth halo background and any structure, we will 
use a Maxwellian distribution, characterized by an rms velocity dispersion <r v , to describe 
the WIMP speeds, and we will allow for the distribution to be truncated at some escape 
velocity v esc , 

( -^—) 3/2 p- 3v2 / 2CT ? for Ivl < v 
n csc {2^ v ) e > lor|v|<t W (15) 

0, otherwise. 

Here 

iV esc = erf (2) - 2zexp(-^)/7T 1 / 2 , (16) 
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with z 



z /vq, is a normalization factor. The most probable speed, 



(17) 



is used to generate unitless parameters such as z. For distributions without an escape 
velocity (t> esc -> oo), N csc = 1. 

The WIMP component (halo or stream) often exhibits a bulk motion relative to us, so 
that 

/(u) = /(v obB + u), (18) 

where v Q b s is the motion of the observer relative to the rest frame of the WIMP component 
described by Eqn. (fl5|) ; this motion will be discussed below. For such a velocity distribution, 
the mean inverse speed, Eqn. (jSJ), becomes 



voy 



i 



2N esc v y 
1 

2N csc v y 

o, 



erf(x+y) — erf(x — y) — A^yz 
erf(z) — erf(x-y) — -^(y + z — x)e 



for z < y, x < \y — z\ 
for z > y, x < \y — z\ 
for \y — z\ < x < y + z 
for y+z < x 



(19) 



where 
recall v 



(20) 



x = v min / v Q , y = v ohs / v , and z = v esc / v ; 

is given by Eqn. (151) . Here, we use the common notational convention of rep- 
resenting 3- vectors in bold and the magnitude of a vector in the non-bold equivalent, e.g. 

Vohs = | V bs | • 

Due to the motion of the Earth around the Sun, v Q b s is time dependent: v t> s = v b s (£). 
We write this in terms of the Earth's velocity V e relative to the Sun as 



VobsW = v + V© [e 1 cosuj(t - ti) + e 2 smuj(t - tj 



(21) 



where u = 27r/year, v is the Sun's motion relative to the WIMP component's rest frame, 
V© = 29.8 km/s is the Earth's orbital speed, and £\ and s 2 are the directions of the Earth's 
velocity at times t\ and t\ + 0.25 years, respectively. With this form, we have neglected the 
ellipticity of the Earth's orbit, although the ellipticity is small and, if accounted for, would 
give only negligible changes in the results of this paper (see Refs. 37|, |51( for more detailed 
expressions). For clarity, we have used explicit velocity vectors rather than the position 
vectors ei and e 2 used in Refs. (HoL 52 ] and elsewhere (the position vectors are more easily 
generalized to an elliptical orbit); the two bases are related by S\ = — e 2 and £ 2 = 

In Galactic coordinates, where x is the direction to the Galactic Center, y the direction 
of disk rotation, and z the North Galactic Pole, 



s 1 = (0.9931,0.1170,-0.01032), 
e 2 = (-0.0670,0.4927,-0.8676) 



(22) 
(23) 



where we have taken e± and £ 2 to be the direction of the Earth's motion at the Spring 
equinox (March 21, or t±) and Summer solstice (June 21), respectively. 

In this paper, we will take the halo to be a simple non-rotating isothermal sphere 0], 
also referred to as the Standard Halo Model (SHM); more complicated halo models will be 
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examined in a future paper. In this model, typical parameters of the Maxwellian distribution 
for our location in the Milky Way are ushm = 270 km/s and v esc = 650 km/s, the latter being 
the speed necessary to escape the Milky Way (WIMPs with speeds in excess of this would 
have escaped the galaxy, hence the truncation of the distribution in Eqn. (115]) ). Unlike the 
Galactic disk (along with the Sun), the halo has essentially no rotation; the motion of the 
Sun relative to this stationary halo is 

V©,SHM = V L SR + v 0jPec , (24) 

where v L sr = (0,220,0) km/s is the motion of the Local Standard of Rest and v 0jPec = 
(10, 13, 7) km/s is the Sun's peculiar velocity. The Earth's speed relative to the halo, t>obs(^)> 
is maximized around June 1. The local dark matter density po is taken to be the estimated 
average density in the local neighborhood, 0.3 GeV/cm 3 . 



C. Annual Modulation 

It is well known that the count rate in WIMP detectors will experience an annual modu- 
lation as a result of the motion of the Earth around the Sun described above 0, 0]. In some 
cases, but not all, the count rate (Eqn. ([1])) has an approximate time dependence 

JD 

— (E, t) « S (E) + S m (E) cosu(t - t c ), (25) 

where t c is the time of year at which v b s {t) is at its maximum. Sq(E) is the average 
differential recoil rate over a year and S m (E) is referred to as the modulation amplitude 
(which may, in fact, be negative). The above equation is a reasonable approximation for the 
SHM we are considering in this paper, but is not valid for all halo models, particularly at 
some recoil energies for dark matter streams; see Ref. |53j for a discussion. For the SHM, 

Sm(E) = \ 

Experiments such as DAMA will often give the average amplitude over some interval, 

S m = - 1 - r dES m (E). (27) 
&2 - -cq Jei 



j|(£, June 1) - ^(E, Dec 1) 



(26) 



D. Parameter Space 

Many of the parameters that factor into the expected recoil rates for a scattering detector 
are unknown, including the WIMP mass, four WIMP-nucleon couplings (SI and SD cou- 
plings to each of protons and neutrons), the local WIMP density, and the WIMP velocity 
distribution in the halo. In this paper, we shall fix the halo model to the SHM and the local 
density to 0.3 GeV/cm 3 . In addition, we shall take f p = f n (equal SI couplings) so that 
there are only three independent scattering couplings; the SI coupling will be given in terms 
of the SI scattering cross-section off the proton, cr p SI . The parameter space we examine will 
then consist of the four parameters m, cr p SI , a p , and a n . 
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Experiment Element Exposure Energies Quenching Efficiency Constraint Ref. 







[kg-day] 


[keVee] 


Factor (Q) 






(t) 




DAMA 


Nal 


2.99 x 10 5 


2-20 


0.3, 0.09 (i) 


1 




(various) 




CDMS 


Si 


6.58 


5-55 


1(§) 


Eqn. 


USD 


Poisson 


f6] 




Ge 


397.8 


10-100 


1(§) 


Eqn. 


J2HD 


Poisson 


m 


CoGeNT 


Ge 


8.38 


0.23-4.1 


0.2 


Eqn. 


JM} 


BP 


[8,9] 


CRESST I 


A1 2 3 


1.51 


0.6-20 


1(§) 


1 




BP 


[12] 


TEXONO 


Ge 


0.338 


0.2-8 


0.2 


Eqn. 


(32) 


MG 


[15,16] 


XENON10 


Xe 


316.4 


6.1-36.5 


1 (§) 


Eqn. 


HMD 


MG 


[11,18] 



TABLE I: Experiments used in this study. Notes to the table: (f) MG is the maximum gap 



method [55J], BP is a Poisson statistics based constraint for binned data as described in the text; 
(t) some portion of recoils may undergo channeling with Q 1 (see text); (§) energies scaled to 
recoil energies. 



III. NULL EXPERIMENTS 



Numerous experiments have searched for a dark matter signal, but all these experimental 
results, apart from DAMA, are consistent with no WIMP signal. Here, we shall examine 
the constraints from these null results; the positive DAMA signal will be addressed in the 
next section. We include here only a few experiments that provide some of the strongest 
constraints at various points in the WIMP mass-coupling parameter space. We further 
restrict ourselves to experiments that provide sufficient data for us to independently generate 
constraints. Notably, this excludes COUPP [l(J 11] for which insufficient information on 
the data runs is provided to both reproduce their results and extend them to arbitrary 
couplings. COUPP is likely to provide strong constraints in areas of parameter space that 
other direct detection experiments do not probe and could significantly affect our results. 
The experiments we include are listed in Table [H 

Most of the experiments observe some recoils; these recoils may be due to background 
events (e.g. neutrons from cosmic ray showers), WIMP scatters, or a combination of the 
two. While tighter constraints may be found by modelling the background contribution, we 
produce conservative constraints by analyzing the experimental results without background 
subtraction. Additionally, some of the data we use (as presented in the literature) is binned, 
whereas analysis of the unbinned data (as available to the experimental groups) generally 
allows for better constraints on the parameter space. Some experiments are also able to 
combine and analyze their data sets more effectively than we can. Where other ambiguities 
arise, we err on the side of caution. For these reasons, our constraints are generally more 
conservative (by as much as a factor of ~3 in some cases) than those provided by the exper- 
imental groups themselves. Readers should refer to the original experimental publications 
for more thorough analyses and, in some cases, the inclusion of background subtraction. 

We will adopt the standard convention of displaying 90% CL. exclusion regions for these 
null experiments. For comparison, however, we will also indicate the parameter space ex- 
cluded at 3a and 5a for a few of the experiments. We will first give a description of the null 
experiments used in our analysis, then we will describe the various statistical methods used 
to provide the constraints for these experiments. 
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A. Descriptions 



CDMS. CDMS II has a large exposure for its Germanium detectors, 397.8 kg-days, with a 
threshold of 10 keV 0]. Observed energies are calibrated to recoil energies, so no quenching 
factor is required. The efficiency of observing nuclear recoils after data cuts is: 

fo.25 + 0.05 (E :l ^ eV) for lOkeV < E < 15keV, 
e(E) = < 5keV (28) 

v ; 10.30 for 20 keV < E < 100 keV. v ; 



While newer data sets exist, we also include the Silicon data set from CDMS-SUF 
due to the low 5 keV threshold, which is lower than the 7 & 10 keV thresholds of the newer 
CDMS II Si data sets. This CDMS-SUF set involves 65.8 kg-days exposure on a 0.100 kg 
Si detector. With an 80% muon anti-coincident cut efficiency, 95% nuclear recoil band cut 
efficiency, and additional energy-dependent cuts, the total efficiency can be approximated 
by 1 : 

e(E) = 0.80 x 0.95 x J" + ^I^E, fa ^ < E < ^ (29) 
V ; [0.40 + 0.10 ( s ° cV V) for 20 keV < E < 100 keV. V ; 

Over a 5-100 keV recoil energy range, two events were detected at 55 & 95 keV that are 
consistent with expected backgrounds. In our analysis, we will use zero observed events over 
a 5-55 keV recoil range. While it is not a statistically sound practice to a postiori select the 
bin to avoid observed events, we are primarily interested in this data set to examine low mass 
( < 30 GeV) WIMPs, for which the observed events are at energies too high to be compatible 
with WIMP scatters. As such, using the chosen bin will result in constraints comparable 
to more formal statistical methods, such as S. Yellin's maximum gap method [55j, at low 
WIMP masses. For higher mass WIMP constraints based solely on this Silicon data, the 
two observed events should be included over e.g. a 5-100 keV range. With inclusion of the 
Ge data, however, the contribution of Silicon to the constraints at higher WIMP masses is 
negligible. 

For the combined CDMS Si+Ge data set, we use Poisson statistics to provide a constraint 
based upon zero observed events. 

CoGeNT. This experiment has a low exposure, 8.38 kg-days, on Germanium detectors, 
but has a very low threshold of 0.23 keVee 0, 0]. Data is binned in observed energies. The 
CoGeNT collaboration has provided us with data in smaller bins than given in their paper; 
the lowest energy bins and the number of observed events in those bins are presented in 
Table HT1 [oj]. From the number of events, exposure, and t§ values in each bin, we estimate 
the efficiency to be: 

^ = °' 66 - 50W' (30) 
where the efficiency is given as a function of observed energy (E' = QE). We take the 
quenching factor to be Qce — 0.20. The widths of the bins are smaller than the energy 
resolution of the CoGeNT detector, so the latter must be taken into account when analyzing 
their results. The energy resolution used by CoGeNT is: 



a(E') = v/(69.7eV) 2 + (0.98 eV)£'. (31) 



1 See also Fig. (3) of Ref. 54|, but note that figure imposes a cut-off at 10 keV and an additional 75% cut 
over 10-20 keV due to problems in one of the Ge detectors; such cuts do not apply in this case. 
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Energy Events Energy Events 

[keVee] [keVee] 



0.338 - 


0.371 


4430 


0.661 - 


0.693 


8 


0.371 - 


0.403 


1104 


0.693 - 


0.725 


5 


0.403 - 


0.435 


140 


0.725 - 


0.757 


18 


0.435 - 


0.467 


29 


0.757 - 


0.790 


15 


0.467 - 


0.499 


12 


0.790 - 


0.822 


12 


0.499 - 


0.532 


12 


0.822 - 


0.854 


5 


0.532 - 


0.564 


19 


0.854 - 


0.886 


12 


0.564 - 


0.596 


9 


0.886 - 


0.919 


16 


0.596 - 


0.628 


12 


0.919 - 


0.951 


9 


0.628 - 


0.661 


8 


0.951 - 


0.983 


8 



TABLE II: CoGeNT binned data. Only the lowest energy bins are included here; an additional 
98 bins over 0.983-4.14 keVee are not shown. Contraints at low WIMP masses mainly come from 
these lower energy bins. 



There are significant background sources that contribute to the large number of observed 
events. The CoGeNT collaboration assumes specific forms for this background and subtracts 
it during their analysis. Here, we make no assumptions on the background and perform no 
background analysis, so our constraint is more conservative. 

For the binned CoGeNT data without background subtraction, we generate constraints 
using the binned Poisson method described below. 

CRESST I. CRESST I used AI2O3 detectors with an exposure of 1.51 kg-days over an 
energy of 0.6-20 keV [12l ]. The energy scale is calibrated to the recoil energy. Only scattering 
off of Aluminum is considered for the spin-dependent scattering case (Oxygen is ignored). 
The energy resolution used by CRESST I is: 

a(E) = v/(0.220keV) 2 + (0.017£) 2 (32) 

Data is binned with a significant presence of backgrounds. We generate the CRESST I 
constraints using the binned Poisson method. 

TEXONO. Another low exposure, low threshold Germanium experiment similar to Co- 



GeNT, TEXONO has an exposure of 0.338 kg-days over energies of 0.2-0.8 keVee |l5l . [16 



The quenching factor is taken to be Qce — 0.2. The combined efficiency after the ACV 



(98.3%), CRV (91.5%), and PSD cuts (energy dependent, see Figure 3 of Ref. [15J]) is esti- 
mated to be: 

e(E>) = 0.983 x 0.915 x {""S^ *» O.lkeV < E' < 0.35 k eV, 
( J [1 for 0.35 keV < E . 

The constraints are generated using the unbinned maximum gap method |55j, which ex- 
amines the likelihood of observing the gaps in energy between observed events. TEXONO 
provides two intervals between event energies, 0.198-0.241 keVee and 1.39-1.87 keVee, that 
are used for the gaps in this method; other gaps in the full 0.198-8 keVee range are ig- 
nored (ignoring other intervals can lead to only more conservative constraints in this type 
of analysis). 
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Some concerns have been raised regarding the claimed behavior of TEXONO's detector 



response at very low energies that calls into question their low energy results [56J. We shall 
not address these concerns here. 

XENON. XENON10 has a large exposure, 316.4 kg-days, on Xenon over 6.1-36.5 keV 



171 . Il8l . l57j |: energies are calibrated to recoil energies (Q = 1). That calibration depends 



upon a quantity C e g, described in their papers, that was initially estimated to be 0.19 



with large uncertainties. More recent results indicate it might be lower [57J and we use a 
conservative value of 0.14. Use of this lower value requires all the nuclear recoil energy scales 
in the XENON10 data to be multiplied by a factor of 0.19/0.14, e.g. the 6.1-36.5 keV energy 
range is determined by multiplying the original range of 4.6-26.9 keV by this factor. We 
approximate the energy dependent efficiency by: 

e(E) = 0.46(1 — ) (34) 

v ' v 135keV ; V ' 

and use an energy resolution of j58[ 



a(E) = (0.579 keV) ^E/keV + 0.02LE. (35) 

XENON10 observed 10 events, with an average background of 7 events expected. We use 
the maximum gap method to generate XENON10 constraints. 

Super-Kamionkande. Aside from the direct detection experiments listed above, we shall 



also show SuperK constraints from WIMP annihilation in the Sun [29j. SuperK determines 
these constraints by searching for high energy neutrino products of the WIMP annihilation 
process at the center of the Sun, where the annihilation rate is dependent upon the capture 
rate of WIMPs and the capture rate depends upon the scattering cross-section for WIMPs 
off of nuclei in the Sun. Due to the large proportion of Hydrogen, the neutrino flux is 
particularly dependent upon the scattering cross-section off of the proton (both SI and 
SD). While there is a weak dependence on WIMP-neutron scattering cross-section through 
the trace amounts of heavier elements, we consider here only constraints for WIMP-proton 
scattering and only in the SD case. The SuperK collaboration also provides SI constraints, 
but only for combined Earth and Solar WIMP annihilations, which we do not include here. 
Their SI constraints do not affect our results. The SD constraint is given at a 90% CL. 
We have not performed our own reanalysis of the SuperK data and simply present the 



constraints as determined by the SuperK collaboration's own analysis (see Ref. [7l| for a 
reanalysis) . 

The SuperK constraints rely on the following assumptions: (1) the WIMP/anti-WIMP 
abundance is not highly asymmetric, which would suppress the annihilation rate; (2) the 
annihilation cross-section is sufficiently high so that the capture rate of WIMPs in the Sun 
(via scattering off of nuclei in the Sun) is in equilibrium with the annihilation rate; and 
(3) the WIMP does not annihilate predominantly into the light quarks, which do not yield 
neutrinos in sufficient quantities and energies to be observed. For the Constrained Minimal 
Supersymmetric Standard Model, these assumptions are mainly satisfied in the parameter 
space of interest. In general, however, these assumptions need not be satisfied, so the SuperK 
constraints should be applied with caution. 

Constraints comparison. Constraints for these null experiments are shown in Figure [T] 
for the case of SI scattering. We show here and in later figures regions excluded at the 
90% level; for comparison, exclusion regions at 3a and 5a are also shown here for a few 
of the null experiments. The regions above the curves are excluded. The high exposure 
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FIG. 1: Contraints on spin-independent (SI) scattering cross-sections for various experiments with 
null results. Cross-sections below each line are excluded by the given experiment at the 90% (solid), 
3a (dashed), and 5a (dotted) confidence levels. The same coloring scheme will be used in later 
figures. 



experiments, CDMS and XENON10, provide very strong constraints at high WIMP masses, 
but provide no constraints at low WIMP masses. Only the low threshold experiments — 
CoGeNT, CRESST I, and TEXONO— provide constraints at low WIMP masses (due to the 
fact that light WIMPs only yield low energy collisions, which high threshold experiments 
cannot detect). 



B. Constraints 

For the various null experiments, we define constraints in parameter space at a certain 
exclusion level 1 — a, typically 90%, as the parameters for which the probability of seeing 
the experimental result is a. Parameters outside the corresponding contours would yield the 
observed result with a probability less than a. We say the parameters within those contours 
are compatible within the 1 — a level and parameters outside are excluded at the 1 — a level. 
The probabilities are determined use different statistical values, described here, in different 
cases. 

Poisson. For data in a single bin, we use Poisson statistics to exclude parameter space 
as follows: parameters that predict an average number of events \i in that bin are excluded 
at a level of 1 — a if the probability of seeing as few as the observed number of events (zero 
in the case of CDMS Si+Ge) is less than a. For zero observed events, that corresponds to 
upper limits on // of 2.3 at a 90% exclusion level, 5.9 at 3a, and 14.4 at 5cr. Here, we use 
Na to indicate the probabilities associated with the ±Na bounds of a normal distribution. 

Binned Poisson (BP). For binned data with unknown background, we exclude parameter 
space using what we shall refer to as a "binned Poisson" (BP) technique, denned as follows. 
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If ctbin is the probability of seeing below a certain number of events in a bin (the lower tail 
of the Poisson distribution), the probability a of seeing at least one bin measurement falling 
below the «bin lower tail of that bin is related by: 

1 _ a = (i _ a bin )^ to) ( 36 ) 

where iVbi n is the number of bins. For a desired exclusion level of 1 — a, we then require at 
least one bin to be at the bin's exclusion level of 1 — abim where «bm can be determined from 
a using the above equation. For example, for a desired exclusion level of 90% in two bins, 
the observed number of events in at least one bin should be excluded at the 94.9% level. In 
other words, if there is only a 5.1% chance of seeing below a certain number of events in 
each bin, there is a 10% chance of the observed number of events in at least one those bins 
to be below that number of events. In practice, the 1 — a exclusion contour is determined 
by, at each WIMP mass, increasing the scattering cross-section until the number of events 
Nk observed in at least one bin is in the 1 — ctbm lower tail of the Poisson distribution 
with average fi^, where fi^ is the theoretical predicted average number of events in that bin. 
Background events make the observed number of events (signal + background) less likely 
to fall into the lower tail of a Poisson distribution based upon the expected signal alone; 
bins with larger signal-to-noise will likely be the first to reach the 1 — abm exclusion level 
(i.e. too few events) and define the overall constraint. Then this technique determines the 
overall constraint from the bin that would provide the strongest constraint, but takes the 
required statistical penalty for choosing that bin, since that bin can simply be one in which 
a statistical fluctuation produces a lower than expected number of events. See Ref. for 
further discussions. This type of analysis is not optimal, but is simple and straightforward to 
perform and works reasonably well for the large background cases of CoGeNT and CRESST I 
here. 

Maximum Gap (MG). For unbinned data with unknown background, constraints are 
generated using S. Yellin's maximum gap (MG) method [55j, which examines the likelihood 
of observing the gaps in energy between observed events. This ubinned method generally 
provides a stronger constraint that the previous (binned) methods for the case of an unknown 
background. 

IV. DAMA 

The DAMA/Nal experiment [2J used a large Nal detector to search for the small modula- 
tion in the WIMP scattering rate and became the first direct detection experiment to observe 
a positive signal. That signal has since been confirmed by the successor to DAMA/Nal, 
DAMA/LIBRA [3] (also a Nal detector). These two experiments remain the only direct 
detection experiments to observe a signal. 

The original DAMA/Nal experiment released data for only two independent bins, a 
low energy bin at 2-4, 2-5, or 2-6 keVee (these are not independent) with a non-trivial 
modulation amplitude and a high energy bin at 6-14 keVee with an amplitude consistent 
with zero. The DAMA/LIBRA results, in combination with the earlier DAMA/Nal results, 
have been released for 36 bins over a range of 2-20 keVee (we shall henceforth use "DAMA" 
to refer to this combined data set). The modulation amplitudes for these bins, taken from 
Figure 9 of Ref. [H, are given in Table IHH The la bounds on the amplitudes for these 
bins are shown as pink boxes in Figure [2j For comparison, modulation amplitudes are 
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Energy 


Average S m 


Energy 




Average S m 


[keVee] 


[cpd/kg/keV] 


[keVee] 




[cpd/kg/keV] 


2.0 


- 2.5 


0.0162 ± 0.0048 


11.0 - 


11. 


.5 


-0.0021 ± 0.0039 


2.5 


- 3.0 


0.0287 ± 0.0054 


11.5 - 


12. 


.0 


-0.0010 ± 0.0040 


3.0 


- 3.5 


0.0250 ± 0.0055 


12.0 - 


12. 


.5 


0.0015 ± 0.0040 


3.5 


- 4.0 


0.0140 ± 0.0050 


12.5 - 


13. 


.0 


0.0020 ± 0.0040 


4.0 


- 4.5 


0.0101 ± 0.0045 


13.0 - 


13. 


.5 


-0.0019 ± 0.0040 


4.5 


- 5.0 


0.0118 ± 0.0041 


13.5 - 


14. 


.0 


0.0015 ± 0.0040 


5.0 


- 5.5 


0.0039 ± 0.0041 


14.0 - 


14. 


.5 


-0.0022 ± 0.0040 


5.5 


- 6.0 


0.0030 ± 0.0039 


14.5 - 


15. 


.0 


-0.0008 ± 0.0040 


6.0 


- 6.5 


0.0059 ± 0.0038 


15.0 - 


15. 


.5 


-0.0021 ± 0.0039 


6.5 


- 7.0 


0.0011 ± 0.0036 


15.5 - 


16. 


.0 


-0.0022 ± 0.0038 


7.0 


- 7.5 


-0.0001 ± 0.0036 


16.0 - 


16. 


.5 


0.0053 ± 0.0037 


7.5 


- 8.0 


0.0004 ± 0.0036 


16.5 - 


17. 


.0 


0.0060 ± 0.0037 


8.0 


- 8.5 


-0.0014 ± 0.0037 


17.0 - 


17. 


.5 


0.0007 ± 0.0037 


8.5 


- 9.0 


0.0039 ± 0.0037 


17.5 - 


18. 


.0 


0.0041 ± 0.0035 


9.0 


- 9.5 


-0.0034 ± 0.0037 


18.0 - 


18. 


.5 


0.0011 ± 0.0035 


9.5 • 


- 10.0 


-0.0071 ± 0.0038 


18.5 - 


19. 


.0 


-0.0014 ± 0.0035 


10.0 


- 10.5 


0.0086 ± 0.0038 


19.0 - 


19. 


.5 


0.0006 ± 0.0034 


10.5 


- 11.0 


-0.0018 ± 0.0039 


19.5 - 


20. 


.0 


0.0056 ± 0.0034 



(f) 10.0 - 20.0 0.0011 ± 0.0008 



TABLE III: DAMA/Nal + DAMA/LIBRA (DAMA) modulation amplitudes for the given observed 
energy bins. Data is taken from Figure 9 of Ref. [3] . (f) The bins over 10-20 keV have been combined 
into a single bin for use with some of the analysis techniques discussed in the text. 

shown in blue boxes for a 2 bin data set: 0.0215 ± 0.0026 cpd/kg/keVee over 2-4 keVee and 
0.0005 ± 0.0010 over 6-14 keVee. 

We use the quenching factors Qno, = 0.3 and Qj = 0.09 in our calculations, with some 
exceptions as discussed in the following section. 

A. Detector Effects 

Several physical and instrumental effects add complications to the basic description of 
experimental recoil rates as described in Section [Til These include a finite detector energy 
resolution and two effects, ion channeling and the Migdal effect, that affect how much of the 
deposited energy goes into the observed channels. 

1. Energy Resolution 

The DAMA/LIBRA apparatus, as with all experiments, cannot determine the event ener- 
gies with perfect precision. If the expected amount of energy in the channel DAMA/LIBRA 
observes is E' = QE, the measured energy will be normally distributed about this energy 
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FIG. 2: DAMA modulation amplitude data as well as best fit spectra for various cases. Boxes 
represent the la bounds of the measured modulation amplitude for each bin: pink for the full 36 
bin data set and blue for a 2 bin data set (2-4 & 6-14 keVee). Also shown are the S m spectra 
for the best fit SI cross-sections (in the SI only coupling case) at a WIMP mass of 4 GeV for the 
2 bin data (blue) and 36 bin data (red), with (dashed) and without (solid) inclusion of the ion 
channeling (IC) effect. The x 2 values (over the degrees of freedom) are indicated in the legend. 
For comparison, the best (SI only) fit at any mass, which occurs around 80 GeV, is shown (green). 



with a standard deviation 60], 



a(E') = (0.448 keV) ^jE'/keV + 0.0091E'. (37) 

While this equation technically applies only to the DAMA/LIBRA apparatus, we use the 
same energy resolution for the combined data of DAMA/Nal and DAMA/LIBRA. We in- 
clude this finite energy resolution in all our DAMA calculations. 



2. Ion Channeling 

Typically, only a small fraction Q of the recoil energy, 0.3 for Na and 0.09 for I, goes into 
a mode that DAMA measures (scintillation). The rest is converted into, e.g., phonons/heat 
as the recoiling nucleus collides with other nuclei and is not observed. As pointed out by 



the DAMA Collaboration 61J] , nuclei that recoil along the characteristic axes or planes of 
the crystal structure may travel large distances without colliding with other nuclei. For 
ions such as recoiling Na + or I~, the recoil energy in this case is primarily transferred to 
electrons (which is observable) rather than other nuclei (which becomes heat). Some recoils 
that undergo this ion channeling (IC) have Q ~ 1. 

DAMA has used simulations to determine the energy-dependent fractions of recoils off of 
Na (fNa) and I (//) that will have Q ~ 1 due to this IC effect. These fractions, given by 
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Figure 4 of Ref . 61[ , are well approximated by 



f Na = (38) 

and 

fj = io-V^/("-Btov) ) (39) 

where these fractions are given as functions of the recozZ energy. 

We shall examine the DAMA results both with and without including the IC effect and 
shall use "IC" in the figures to indicate those analyses that include it. We have not included 
this effect in other experiments where it may also be applicable, such as CoGeNT and 
TEXONO 2 . 



3. Migdal Effect 

In the 1940's, Migdal pointed out that the rapid change in velocity of a recoiling nucleus 
can cause bound electrons to become excited or ionized 62|. In the DAMA detector, the 



energy of these ionized electrons goes into observable modes. At the same time, the apparent 
recoil energy is reduced as the energy is lost to the ionized electron almost immediately, 
before the nucleus travels far enough to interact with other nuclei and electrons. This 
Migdal effect (ME) can then affect the amount of observed energy coming from a scattering 
event at a given energy. Due to the complicated nature of this effect, we shall not show its 
implementation here, but it results in a rather complicated response function $ that may 
be used in Eqn. (J5]). The implementation of the Migdal effect in an analysis of DAMA is 
fully described in Ref. 63[. We do not include the Migdal effect in our analysis at this time, 



but we plan to include it in a future paper. 



B. Analysis Techniques 

There are many possibilities for analyzing the DAMA results; the choice of statistical 
tests depends upon how one wishes to interpret the results. Here, we will use three different 
tests to produce different regions in parameter space: a likelihood ratio method with a global 
fit of four parameters (which we refer to as the "likelhood ratio" method from now on) to 
find the preferred (best fit) parameters to produce the DAMA signal, a likelihood ratio 
method with a one parameter fit in raster scans to find the preferred scattering cross-section 
to produce the DAMA signal at each WIMP mass (which we call the "raster scan"), and a 
X 2 goodness-of-fit test (g.o.f.) to indicate which parameters are compatible with the DAMA 



signal. See Ref. 64] for a short review of statistics or Ref. [65] for more extensive discussions. 



2 Channeled events in experiments that measure phonons/heat, such as CDMS, would not be seen as 
any events with relatively high ionization/scintillation and low heat are characteristic of electromagnetic 
backgrounds (e.g. 7-rays), which induce electron rather than nuclear recoils, and are thrown out. 
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Mass [GeVJ 


r ui 
•^Sl [pbj 


a P (o-p,SD [pb]) 


a n (<7n,SD [pbj) 


Xmin ( d -°' f -) 


Global 


81.0 


2.0 x 10 


0.0 


(0.0) 


0.0 (0.0) 


27.1 


(32) 


w/ channeling 


78.6 


2.7 x 10~ 5 


0.62 


(0.13) 


0.70 (0.17) 


25.4 


(32) 


hi only 


81.0 


2.0 x 10~ 5 





(0) 


(0) 


27.1 


(34) 


/ l 1 ■ 

w/ channeling 


80.0 


3.0 x 10~ 5 





(0) 


(0) 


26.3 


(34) 


SD only 


9.1 





14.2 


(59.4) 


-176 (9020) 


29.3 


(33) 


/ l 1 ■ 

w/ channeling 


12.0 





1.07 


(0.35) 


-12.8 (50.2) 


30.9 


(33) 


bD proton only 


11.5 





1.70 


(0.88) 


(0) 


30.8 


(34) 


/ l 1 ■ 

w/ channeling 


11.4 





1.33 


(0.53) 


(0) 


31.0 


(34) 


OTA j_ 1 

bD neutron only 


11.0 








(0) 


20.3 (125) 


31.0 


(34) 


w / cnaiiiieiiiig 


11. y 


n 





(0) 


7 8 CIS fil 


30.9 


(34) 


SD only (a n = a p ) 


11.5 





1.57 


(0.75) 


1.57 (0.75) 


30.8 


(34) 


w/ channeling 


11.4 





1.15 


(0.40) 


1.15 (0.40) 


31.0 


(34) 


SD only (a n = — a p ) 


11.6 





1.88 


(1.05) 


-1.88 (1.05) 


30.8 


(34) 


w/ channeling 


11.2 





1.56 


(0.73) 


-1.56 (0.73) 


31.1 


(34) 



TABLE IV: Parameters at which the DAMA x 2 is minimized. The global minimum is found by 
allowing all parameters (m, Cp,sij a p, and a Q ) to vary, while the other cases allow only a subset 
of these parameters to vary freely with the other parameters fixed to zero or fixed in terms of the 
free parameters. Bold parameters in the table are those that have been fixed. 



1. Likelihood Ratio 

To determine the most likely parameters for producing the DAMA signal, we use the 
maximum likelihood method, based on the likelihood ratio 

L(S mt k\m, cr P;SI , Op, q n ) 

r/rr I - * * * \ 5 V^UJ 

where L is the likelihood function, S m ^ are the observed modulation amplitudes in each bin, 
and m, <x p ,si) a p , and a n are the values of the parameters that maximize the likelihood for 
the observed S m ^- The denominator of the above equation is the maximum likelihood value 
L max . Confidence regions in the parameters are determined by 

2 In L(m, er P)SI , a v ,a n ) > 2 lnL max — 2 A InL, (41) 

where the value of A In L corresponds to the confidence level (C.L.) of the confidence region. 
Since each bin is normally distributed, this equation may instead by given in terms of the 

x\ 

X 2 (m, a PjS i, Op, On) < xlain + A X 2 , (42) 

where 

( q _ cTh \ 2 

X {™>, <7 p ,si, Op, On) = 2 ' ( 43 ) 

k G k 

(Tfc is the uncertainty in S m ^ (see Table [TTT]) . Sj r ^ k = S^ k {m, c p ,si, a p , a n ) is the expected 
amplitude in a particular bin for the given parameters, and Xmin i s the minimum value of 
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X 2 ■ The Xmin an d the parameters at which it occurs is given in Table HVl in this global (four 
parameter) case both with and without IC. Also given are the minima in various specific 
cases, such as SI only scattering. 

The Ax 2 limit corresponding to different C.L.'s depends upon the number of parameters 
that have been minimized and, in the large data sample limit, is determined from a x 2 dis- 
tribution with the degrees of freedom (d.o.f.) equal to the number of minimized parameters. 
For the four parameters considered here, Ax 2 = 7.8 (90% C.L.), 16.3 (3a), 34.6 (5a), and 
60.3 (7a). 

The confidence region for DAMA as described here is, in fact, a 4-dimensional region in 
the (m, a P) si, a p , a n ) parameter space. We shall be showing 2-dimensional slices of this larger 
region for particular cases, such as the a p = a n = slice corresponding to SI only scattering. 
This is not equivalent to fixing a v = a Q = and determining a 2-dimensional confidence 
region by minimizing only over (m, a Pi si). In the latter case, a confidence region will always 
be defined in a 2-dimensional plane while, in the former, there may be no slice in that plane 
if that plane represents a poor "fit" relative to the overall parameter space. 

Note, as will be discussed in Section [V] the confidence regions obtained via this method 
yield the most preferred (best fit) parameters for producing the signal. This method, in 
and of itself, does not imply parameters outside the confidence regions are necessarily a bad 
fit to the data (or, conversely, that parameters inside these regions are a good fit) and one 
should be careful using these regions to compare with other experiments 3 . 



2. Raster Scan 

Since previous papers have examined the possibility of low mass WIMPs as the source 



of the DAMA signal [20J, [21|, [22|, we include an analysis that gives the preferred scattering 
cross-sections at a given WIMP mass. If one specified how to determine the preferred cross- 
sections for an arbitrary WIMP mass, the band formed by these preferred cross-sections 
over all masses is what is referred to as a raster scan. 

The raster scans here involve only two parameters, m and a, where the WIMP mass m is 
the parameter we scan over and a is some scattering cross- sect ion, e.g. the SI cross-section 
(we shall be examining several different cases). At a given m, the preferred values for a 
are determined using the likelihood ratio method just described but leaving one parameter, 
which we call here a, as a variable and fixing the other two couplings/cross-sections to some 
value, usually zero. Thus the x 2 is minimized over a at that particular mass, XRs,min( m ) = 
X 2 {m,d-). The confidence interval in a at m is then 

X 2 (m, a) < XRs.minM + Axrs> (44) 

where limits on A^g, determined from a x 2 distribution with 1 d.o.f. (large data sample 
limit), are 2.7, 9, 16, and 25 at C.L.'s of 90%, 3a, 5a, and 7a, respectively. The band of 
confidence intervals in a over all masses defines a confidence region in m and a at the same 
C.L. as the intervals; we call this region the raster scan region. 

The above formulation is not technically correct when physical boundaries are present, 
such as the requirement of a non-negative cross-section, but it is approximately correct if 



In statistical parlance, the determination of the confidence region for this method is decoupled from the 
goodness-of-fit. 
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the obtained confidence intervals are not near the boundary. As will be seen later, the raster 
scans do yield confidence intervals that cross the physical a = boundary; however, these 
occur only at high WIMP masses that will not be of interest. The raster scans over the 
masses of interest are far enough from the boundary to be approximately correct. Formally, 
one should use a method that can appropriately account for such boundary conditions, such 
as that of Feldman and Cousins [66[] 4 . 




The raster scan itself, by definition, provides no constraint on the WIMP mass and the 
raster scan band includes WIMP masses at which no scattering cross-section is capable of 
reasonably producing the DAMA signal. One may determine which WIMP masses provide 
poor fits to the DAMA data by applying a x 2 goodness-of-fit test at the best fit cross-section 
a at the mass to be examined. The x 2 values at these minima should follow a x 2 distribution 
with 35 d.o.f. (36 bins with a fit to 1 parameter (a)). The x 2 f° r 35 d.o.f. is expected to 
fall below 46.1, 62.8, and 91.7 with probabilities of 90%, 99.73% (3a), and 99.999943% (5a), 
respectively; then the x 2 should exceed those values only 10%, 0.27% and 5.7x 10~ 5 % of the 
time, respectively. 



To conservatively indicate parameters that are compatible with the DAMA data set, we 
also indicate regions at which the x 2 f & U s within a given level using a simple x 2 goodness- 
of-fit (g.o.f.) test on the data. In contrast to the previous two analyses methods, there is no 
fit to the parameters here. The g.o.f. regions are defined as those parameters for which 



where Xgof * s the value at which the x 2 cumulative distribution function (CDF) for e.g. 
36 d.o.f. is equal to the desired level of compatibility. That is, for a desired compatibility 
level of 1 — a, there is a probability a that x 2 wm exceed Xgof- Alternatively, we can say 
parameters outside of the region are excluded at the 1 — a level. 

For the SHM, the expected DAMA modulation amplitudes for energies over 10 keVee 
are negligibly small (statistically equal to zero) over the parameter space examined in this 
paper. The DAMA modulation amplitude data is then expected to randomly vary about 
zero in the higher energy bins, which increases the x 2 an d dilutes the power of the g.o.f. 
test. To improve the ability of the g.o.f. test to exclude some parameters, we combine all 
DAMA bins over 10-20 keVee into a single bin (0.0011 ± 0.0008 cpd/kg/keVee), resulting in 
a total of 17 data bins used for this test. The previous two analyses (likelihood ratio and 
raster scan) are not affected by the numerous high energy bins varying about zero and use 
the full 36 bin data set. For the 17 d.o.f. x 2 distribution, values of 24.8, 37.7, and 61.6 are 
excluded at the 90%, 3a, and 5a levels, respectively. 

While the g.o.f. regions defined here are, in fact, confidence regions — indicating the likely 
parameters to produce the DAMA signal in our theoretical framework — we do not use these 
regions in that manner. The determination of a confidence region assumes that the theoret- 
ical framework is correct, i.e. there exists some choice of parameters that is correct. That 




3. Goodness- of- Fit 



X 2 (m, a PiS i, dp, an) < x 2 gof, 



(45) 



4 The result of using this more formal method will be to shift upwards the upper bounds on a at massses 
where the confidence interval is near the a — boundary, but the shift will be less than a factor of 2. 
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assumption may not be valid if, e.g., the standard halo model is not a reasonable approx- 
imation of the actual halo. Instead, we will more conservatively only use this g.o.f. test to 
exclude parameters outside of the corresponding regions. That is, parameters outside of the 
DAMA g.o.f. regions are incompatible with the DAMA signal at the given level. 

C. Total Rate 

While we mainly examine the DAMA experiment by analyzing the modulation amplitude 
S m , DAMA can additionally constrain parameter space using the total rate (So in Eqn. (1231) ) 
as shown in Figure 1 of Ref . [sj . The DAMA detectors do not strongly discriminate between 
WIMP scatters and background events, so an unknown and possibly large portion of the 
total observed events may be due to backgrounds (the backgrounds are presumed not to vary 
with time and do not contribute to the modulation). In the same manner as the various null 
experiments, we use the total number of events observed by DAMA in 0.25 keVee width bins 
over 2-10 keVee to constrain the parameter space by excluding regions that would predict 
more events than observed. While DAMA shows data for several bins below 2 keVee, we 
conservatively use only the data above their 2 keVee threshold as the behavior of their 
detector at low energies may not be well understood. Due to the large number of observed 
events (in excess of 10,000 in each bin), statistical fluctuations are relatively small and any 
choice of parameters that predicts an So value exceeding the measured rate by more than a 
few percent is incompatible with the data to a very high level. 

As with the null experiments, we will show regions excluded at the 90% level due to 
the total rate for DAMA. Constraints will be show with and without the channeling ef- 
fect. Note that some choices of parameters that produce the correct modulation signal may 
actually predict a total rate larger than observed and, hence, those parameter choices are 
incompatible with the full DAMA data set. 

V. RESULTS 

Here, we apply the statistical techniques described in the previous section to the DAMA 
data and make comparisons with the null experiments. We begin by examining how the 
increase in DAMA data bins affects the results, then we examine the various DAMA detector 
effects, and finally we compare the various experiments. 

A. DAMA Binning 

The DAMA/Nal data was only presented for two independent bins, so our earlier analyses 
(performed prior to DAMA/LIBRA) used only one or two bins. To demostrate how the 
increase in the number of DAMA data bins from 2 to 36 has considerably increased the 
ability to constrain parameter space, we first examine a 2 bin data set. Here, we use the 
2 bins from the combined DAMA/Nal + DAMA/LIBRA data set as given in the previous 
section (2-4 & 6-14 keVee); this may be directly compared with the 36 bin DAMA results, 
which are derived from the same raw data set. The bins and the la bounds on the amplitudes 
are shown as blue boxes in Figure [2j The low energy (2-4 keVee) bin has a positive amplitude, 
while the high energy (6-14 keVee) bin is consistent with no modulation. 
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FIG. 3: Most likely parameters to produce the DAM A signal as determined using 2 bins (2-4 & 
6-14 keVee; blue) and 36 bins (pink) for the case of SI scattering. The dark/light solid regions 
correspond to 90% /3a confidence regions; these regions are surrounded by 5a and 7 a contours. 
The 2 bin data set provides very little constraint on the WIMP mass, while the 36 bin constrains 
the WIMP mass to two regions around 12 GeV and 80 GeV, corresponding to predominantly 
scattering off of Na and I, respectively. Null experiment constraints as described in Figure Q] are 
shown at 90% exclusion levels. 

In the SHM, there are very few recoils in Nal above 6 keVee, so the expected modulation 
amplitude in the high energy bin will nearly always be small. In that case, the upper bin 
provides essentially no constraint on parameter space and there is effectively only one bin 
(the low energy bin) that yields constraints. Since the amplitude is proportional to the 
scattering cross-section, as long as a positive modulation is predicted for the low energy 
bin at some WIMP mass, one can always find a cross-section that yields an excellent fit to 
the data. Only at very low WIMP masses, where the recoils are all of low energy and fall 
below the 2 keVee threshold so there is zero predicted modulation, and at very high WIMP 
masses, where the modulation amplitude is negative so that a good fit requires an unphysical 
negative cross-section, we are unable to fit the 2 bin data. Good fits can be found over a 
broad range of WIMP masses, ~ 2-100 GeV, as indicated in Figure [3] and Figure HI the latter 
of which includes the channeling effect. Effectively, analysis of the 2 bin data requires only 
a fit to the amplitude of the modulation spectrum, not the shape of the spectrum, and only 
weakly constrains the WIMP mass. We used here two data points and four parameters, thus 
the value of the minimum % 2 , which is zero for several points in parameter space, does not 
provide an indication of the goodness of he fit. 

The 36 bin data set breaks the full energy range into much smaller bins; the positive 
signal is now found in six separate bins over 2-5 keVee. To provide a good fit to the data in 
this case, the predicted spectra must not only have the correct amplitude, it must have the 
correct shape to produce the proper signal in all six of those bins. [The higher energy bins 



23 




M WIMP (GeV) 



FIG. 4: Same as Figure El but including the channeling effect for DAMA. The 36 bin low mass 
region around 12 GeV is predominantly due to scatters off of I that undergo the channeling effect, 
as opposed to scatters off of Na as in the no ion channeling case shown in Figure [3j ordinary Na 
scatters and channeled I scatters coincidentally provide good fits at similar WIMP masses. 

are consistent with zero modulation, but again that is typically expected for the spectra in 
the SHM.] Since the shape of the spectrum depends significantly on the WIMP mass, the 
36 bin data set is able to much better constrain the WIMP mass than the 2 bin data set, 
as shown in Figures [3] & HI The reason for the difference can more clearly be seen by the 
spectra in Figure [2J At a low WIMP mass of 4 GeV, the (no channeling) spectra at the best 
fit SI cross-section is shown in blue for 2 bins, with a x 2 of 0.25 on 1 d.o.f., and in red for 
36 bins, with a x 2 °f 87.0 on 35 d.o.f.; spectra including the channeling effect are shown in 
dashed lines with the same colors (with corresponding x 2 given in the legend). At such a 
low WIMP mass, the spectra are very sharp and drop off rapidly as light WIMPs are not 
capable of inducing high energy recoils. The fit to the 2 bin data is, in fact, nearly a perfect 
fit to the low energy bin as the very large amplitude just above the 2 keVee threshold and 
small amplitude over the rest of the 2-4 keVee bin gives the correct average modulation over 
that bin. On the other hand, the fit to the 36 bin data is very poor: it overestimates the 
modulation amplitude in the first bin, but underestimates it for the next five. The shape 
of the spectra required by the 36 bin data cannot be matched at this WIMP mass. For 
comparison, the spectrum is shown in green for the best fit mass, near 80 GeV, which has 
a x 2 °f 27.1 on 34 d.o.f. This spectrum is a very good fit to both the 2 bin and 36 bin data 
sets. 

Of note is the fact that the first (lowest energy) bin actually has a lower amplitude than 
the next few bins. At low enough energies, the modulation actually reverses phase and the 
amplitude becomes negative. The lower first bin may be an indication that the DAMA 
threshold is just above the energy at which the modulation reduces phase. Since that 
energy is dependent upon the WIMP mass, the lowest bin makes an important contribution 
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to constraints on the WIMP mass. 

For 36 bins, there are two mass regions that allow a good fit to the data: around 10- 
15 GeV and 60-100 GeV, as seen in Figure [31 In the lower mass region, scattering off of Na 
produces a spectrum that matches the data, while in the higher mass region, scattering off of 
I produces the appropriate spectrum. In those two regions, scattering is also predominantly 
off of Na and I, respectively. 

The 36 bin data will be used for the remainder of the paper. 

B. DAM A Detector Effects 

Energy resolution. The finite energy resolution serves to smear out the modulation spec- 
trum. It leads to a mild improvement in sensitivity to low mass WIMPs, but is not as 
significant an effect as the two discussed below. The finite energy resolution is included in 
all our results. 

Ion channeling (IC). The channeling effect reduces DAMA's recoil energy threshold. 
Whereas a scattering event with the minimum observed energy of 2 keVee is conventionally 
assumed to be the quenched energy of a 6.7 keV recoil off of Na (Qnu = 0.3) or a 22.2 keV 
recoil off of I (Qi = 0.09), it could also be due to a channeled event with a true recoil energy 
of 2 keV. This lower recoil energy threshold provides greater sensitivity to low mass WIMPs. 
This is apparent when comparing the 2 bin result of Figure H] (channeling) with Figure [3] 
(no channeling), as a much smaller cross-section is sufficient at low masses to produce the 
necessary signal in the channeling case. 

There are four groups of scatters — off of Na and I, each with or without channeling — that 
contribute to the modulation spectrum and each of those groups will have a best fit to the 
data at different regions of parameter space. Whereas there were two "preferred" regions 
without channeling (Figure [3]), there are now potentially four. Predominantly scattering off 
of Na and I, without channeling, should again represent a good fit to the data around 10- 
15 GeV and 60-100 GeV, respectively. The best fit regions to channeled scatters off of Na 
is at even lower masses and is apparent as the bump in the 2 bin region of Figure H] around 
2-4 GeV. However, as evident by the lack of 36 bin regions (only the la contour covers that 
mass) , channeled events off of Na do not actually provide a reasonably shaped spectrum at 
these light masses, even though that is where their "best" spot is. The only other region 
is that due to channeled events off of I. That coincidentally also occurs around 10-15 GeV, 
the same as unchanneled Na events. In fact, the channeled I events provide a much larger 
contribution to the modulation at those WIMP masses than the unchanneled Na events do, 
as evidenced by the reduction in the best fit cross-section by over an order of magnitude from 
Figure [3] to Figure HI Thus, there are two best fit mass regions when channeling is included: 
a ~10-15 GeV low mass region due to predominantly channeled scattering events off of I and 
a ~60-100 GeV high mass region due to predominantly unchanneled scattering events off 
of I. Though the best fit masses without including channeling are nearly the same as when 
including this effect, that is mainly a coincidence: the low mass region in the former case is 
due to scatters off of Na, while the low mass region in the latter case is due to channeled 
scatters off of I. 

Migdal effect. We do not include the Migdal effect in our analysis at this time, but we 
plan to include it in a future paper. 
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FIG. 5: Experimental constraints and DAMA best fit parameters for SI only scattering. The 
DAMA best fit regions are determined using the likelihood ratio method with (green) and without 
(orange) the channeling effect. 



C. DAMA Analysis Techniques 

The likelihood ratio analysis of DAMA yields a 4-dimensional confidence region over 
the (m, er p ,si, %>, o-n) parameter space. In Figures El El and d we show the SI only (a p = 
a n = 0), SD proton-only (cr p SI = a n = 0), and SD neutron-only (cr p SI = a p = 0) slices 
of this confidence region, respectively, with (green) and without (orange) the channeling 
effect. Later, we will show slices in other, mixed coupling cases. This analysis indicates the 
parameters most likely to produce the DAMA signal. If some 2-dimensional plane in the 
full parameter space allows for only a poor fit relative to the best fit point (the global x 2 
minimum), there will be no slice in that plane. The fact that we have found slices in all of 
the SI only, SD proton-only, SD neutron-only, and mixed SD coupling cases is an indication 
that each of these cases is a relatively reasonable possibility for producing the DAMA signal. 
That is, even though we allow for mixed SI and SD couplings, the DAMA signal can still 
reasonably be produced by SI only scattering; the same may be said of the other coupling 
possibilities. One may see this from the minimum \ 2 given in Table IIV1 although the SI 
only case provides a mildly lower Xmin (27.1 without channeling) than that the SD cases 
(29.3-31.0), they are all very similar to the global minimum (27.1). 

By finding a 4-dimensional confidence region, direct comparisons and (to some degree) 
interpolations may be made between the DAMA regions of different coupling cases. Alterna- 
tively, we could have considered each case to be its own 2 parameter model in (m,a), where 
a represents the type of cross-section in that case {e.g. o" Pj si), and found the corresponding 
2-dimensional confidence regions. However, each analysis would be based on a different 
theoretical framework and would not be statistically correlated; one must then be cautious 
when making comparisons between, e.g., the resulting SD proton-only and SD neutron-only 



26 



io 4 



10" 



Oh 



10" - r 
io-Sr 
io- 2 
l o 

10" 




spin-dependent 
(a n = 0, proton only) 



io 1 



10 J 





DAMA (7<x/5o-) 




DAMA (3o-/90%) 




DAMA (7o-/5tr) 
with channeling 
DAMA (3<r/90%) 
with channeling 




CRESST 1 




TFYOMO 

1 JjAUl* V./ 




CoGeNT 




Super-K 




XENON 10 




CDMS I Si 
CDMS 11 Ge 



WIMP 



(GeV) 



FIG. 6: Experimental constraints and DAMA best fit parameters for SD proton-only scattering. 
The DAMA best fit regions are determined using the likelihood ratio method with (green) and 
without (orange) the channeling effect. The CoGeNT and TEXONO constraints are too weak to 
fall within the shown region. 
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FIG. 7: Experimental constraints and DAMA best fit parameters for SD neutron-only scattering. 
The DAMA best fit regions are determined using the likelihood ratio method with (green) and 
without (orange) the channeling effect. Super-Kamiokande (SuperK) provides no constraints in 
this case and is not shown. 
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FIG. 8: Two statistical analyses of DAMA (best fit and raster scan) vs. null results of other 
experiments for SI scattering. DAMA best fit regions (likelihood analysis) are shown by the 
orange regions (90% /3a) and contours (5a /7a). [Technically, these regions are the 2D a p = a n = 
slices of the full 4D confidence region as described in the text.] The blue bands and contours are 
raster scans over the WIMP mass, corresponding to the most likely cross-sections at each WIMP 
mass The data is not used to constrain the WIMP mass in a raster scan analysis; however, the 
brown/orange/yellow line shows the cross-section that minimizes the \ 2 (best fit cross-section) at 
each mass, with the color indicating a goodness-of-fit within the 5a /3a/ '90% level at those cross- 
sections (dotted yellow indicates a goodness-of-fit below 5a, a very poor fit). Null experiment 
constraints are as described in Figure [H The constraint imposed by the total event rate in DAMA 
(as opposed to the modulation) is shown by the dark dotted curve. The likelihood ratio analysis is 
decoupled from a goodness-of-fit test and the lack of an overlap between the DAMA best fit regions 
(as determined by the likelihood ratio) and the allowed regions of other experiments does not, by 
itself, indicate an incompatibility between the experiments. Such comparisons, however, may be 
made using the x 2 goodness-of-fit contours shown in Figure [TUJ The constraint due to the total 
event rate in DAMA (dotted line) shows that some parameters that yield the correct modulation 
amplitude predict too many total events and are therefore excluded. 



confidence regions. In theory, such 2-dimensional confidence regions may be very different 
from the slices of the 4-dimensional confidence region. In practice, for the coupling cases 
considered here, they turn out to be somewhat similar, which is a consequence of the similar 
Xmin i n the different cases. 

As discussed previously, there are generally two preferred WIMP mass ranges: a low 
mass range ~10-15 GeV (due to scatters primarily off of Na or channeled I scatters, when 
channeling is included) and a high mass range ~60-100 GeV (due to scatters primarily off 
of I). Low, high, and intermediate WIMP masses generally produce worse fits (though not 
necessarily bad fits). Previous analyses have found low WIMP masses, below these preferred 
ranges, to be of interest. In addition, the LHC or other (non-direct dectection) experiments, 



28 




10 -8 _| , , , , , , ■ ■ ■ 

10° 10 1 10 2 10 3 



M WIMP (GeV) 



FIG. 9: Same as Figure El but including the channeling effect for DAMA. 

if they observe e.g. neutralinos, might be able to provide the mass of the WIMP. For these 
cases, it is useful to use the raster scan, which does not use the data to make assumptions on 
the WIMP mass and provides appropriate confidence intervals in the scattering cross-section 
for a given (known) WIMP mass. 

The raster scan is shown in Figure [8] (without channeling) and Figure [9] (with channeling) 
for SI only scattering. Regions shown in dark/light blue are the confidence regions at 90%/3<r 
C.L., while the light blue contours indicate the confidence regions at a ba /7a C.L. (the outer 
contour is 7 a). For comparison, the DAMA likelihood ratio regions and null experiment 
constraints are also shown; these regions and constraints are equivalent to the ones shown 
in Figure The raster scan regions pass through the likelihood ratio regions (as should be 
expected). At masses where the raster scan regions/contours extend to the bottom of the 
figure {7a /5a at low masses, 7a at intermediate masses, and for all C.L.'s at high masses), 
the regions actually extend all the way to a zero cross-section, i.e. a cross-section of zero 
falls within the corresponding C.L. at those masses. 

The brown/orange/yellow line along the raster scan indicates the best fit cross-sections 
minimizing the x 2 ) at each WIMP mass. At a particular mass, the reasonableness of the 
best fit cross-section (and corresponding confidence interval) may be tested against the x 2 
distribution with 35 d.o.f. (36 bins with a fit to one parameter): the brown/orange/yellow 
color of the line indicates the best fit cross-section at that mass falls within a 90%/3cr/5cr 
probability level (dotted yellow is below 5<r). Changes in the color often occur near one 
of the likelihood ratio contours; this is only coincidental as the two types of regions use 
different statistical methods and are not in direct correspondence. 

The most appropriate use of the raster scan is if one has an interest in a particular WIMP 
mass (or specific range of WIMP masses) : the raster scan indicates the confidence interval 
(most preferred or best fit values) in the cross-section at that mass and the color of the 
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central line indicates the goodness-of-fit of that particular model/interval 5 . If there is no 
interest in specific WIMP masses, the likelihood ratio method should be used instead. 

The regions determined using the likelihood ratio and raster scan methods indicate the 
best fit parameters for producing the DAMA signal, assuming that signal is actually pro- 
duced by some particular set of parameters in the assumed theoretical framework. These 
analyses, however, are decoupled from the goodness-of-fit: the resulting regions do not in- 
dicate if those parameters are a good fit to the data or a poor one (which may indicate the 
assumed theoretical framework, e. g. the SHM, is incorrect). That is, the best fit parameter 
regions are not equivalent to the regions that provide a good fit to the data. Some parame- 
ters that fall outside of the best fit regions may actually provide reasonably good fits to the 
data (just not as good as parameters within the region), while some parameters inside the 
best fit regions may provide poor fits to the data (just not as bad as parameters outside the 
region). In the likelihood ratio method, only one of these two cases is possible, depending 
upon the value of Xmin' ^- e - there will be either good fits outside the region or bad fits within, 
but not both. 

Because these two regions indicate best fit parameters, but not the "allowed" parameters, 
they should not be used to determine the compatibility with other experiments. For example, 
in the SI case with channeling shown in Figure [5j the XENON10 constraint essentially rules 
out the 4cr DAMA likelihood ratio best fit parameter region (not shown, but falls midway 
between the 3a and 5a regions). However, at about m = 8 GeV and a Pt si = 1CT 5 pb, 
the predicted recoil and modulation spectra are compatible with XENON10 and DAMA at 
about the 90% and 2a levels, respectively, using a \ 2 distribution for the latter. 

The best way of determining the compatibility of the various experiments is to do a 
combined analysis, such as a likelihood ratio analysis, that takes into account all the available 
experimental data 6 . However, this is fairly difficult to do due to the wide range of analyses 
used by the various experiments and the lack of background subtraction makes it even more 
difficult (the likelihood is ill defined in this case). Rather, we will more conservatively 
compare experiments using a \ 2 goodness-of-fit (g.o.f.) test for DAMA for each set of 
parameters (there is no fit to the parameters, so we use a x 2 distribution with 17 d.o.f., 
recalling the high energy bins have been combined for the g.o.f. test). The \ 2 g.o.f. regions are 
shown as black contours in Figure (TU] (without channeling) and Figure [TT] (with channeling) 
where the DAMA data is compatible with the predicted spectrum at the 90%, 3a, and 5a 
level. Parameters outside those contours are incompatible with the data at at least the 
corresponding level; these contours are constraints on the parameter space. These contours 
are appropriate for comparing DAMA with the other experiments since, if a set of parameters 
provides a poor fit to the data of one experiment, including DAMA, it almost certainly will 
be a poor fit in a combined analysis as well. These x 2 g-o-f- contours are a conservative 
means of examining the compatibility as a combined analysis could possibly provide stronger 
constraints, but is unlikely to provide more lenient constraints. 

As they are both determined from the \ 2 °f the data, the contours for the likelihood 



5 A technical point: the raster scan confidence region is the band of confidence intervals over all masses. 
The goodness-of-fit test at the best fit cross-section at a given mass is only a test of the confidence interval 
at that mass. If this test is used to "cut" the band at some WIMP masses, the resulting region would not 
formally be a confidence region. 

6 The best fit parameter region determined via the likelihood ratio method may shift significantly with 
these other experimental constraints added, yet still provide a reasonable fit to the data. 
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FIG. 10: Two statistical analyses of DAMA (best fit and goodness-of-fit (g.o.f.)) vs. null results 
of other expeiments for SI scattering. The best comparison with other experiments is via the 
g.o.f. regions shown here. The most likely parameters to produce the DAMA signal, as determined 
using the likelihood ratio method, are indicated by the orange regions (90%/3<r) and contours 
(5a/7a). The black contours correspond to a conservative x 2 g-o.f. test to the DAMA data. The 
constraint imposed by the total event rate in DAMA is shown by the dotted curve. Null experiment 
constraints are as described in Figure [5j Both DAMA regions shown are based on the x 2 value and 
thus have similar contours shapes, but they have different x 2 cutoffs at a given confidence level and 
are not equivalent; e.g. the 5a likelihood ratio contour is nearly the same as the 3a g.o.f. contour. 
Direct comparisons between DAMA and the other experiments should be done using the x 2 g-o.f. 
regions, not the likelihood ratio regions. 

ratio and g.o.f. methods are similar in shape, as can be seen in Figures [10] & [TlJ However, 
the x 2 corresponding to a given CL. is different in the two methods. In particular, while 
some contours of the likelihood ratio analysis seem to closely match contours of the x 2 g-o.f. 
analysis, they are not of the same level. For example, the 5a likelihood ratio contour is very 
similar to the 3a g.o.f. contour. 

To summarize, the choice of the DAMA analysis method to use should be the following: 
if the most likely parameters to produce the DAMA signal are desired, use the likelihood 
ratio; if interested in a particular WIMP mass or range of masses, use the raster scan; and if 
direct comparisons (particularly, the compatibility) are to be made with other experiments, 
use the x 2 g-o.f. analysis. 

D. DAMA Total Rate 

The total event rate in DAMA imposes constraints on the DAMA compatible parameter 
space; this total rate constraint is shown as a dotted curve in Figure [8] and the following 
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FIG. 11: Same as Figure [T0| but including the channeling effect for DAMA. 

figures. Parts of parameter space that would produce the DAMA modulation signal, such 
as the higher WIMP mass regions for the SI case shown in Figure [10] (without channeling) 
and Figure [TT] (with channeling), predict far more total events in the DAMA detector than 
observed and are therefore ruled out as an explanation for the DAMA signal. 

In general, the total rate provides a strong constraint for the higher WIMP mass region 
that is favorable for the modulation signal in the various coupling cases, excluding the 90% 
(without channeling) and 3<r (with channeling) g.o.f. regions. Thus, high mass WIMPs are 
disfavored as the source of the DAMA signal for the models considered here. The lower 
WIMP mass regions, however, are only minimally constrained by the total rate in DAMA. 

E. Spin-Independent Scattering 

The DAMA best fit parameters and null experiment constraints for SI only scattering are 
shown Figure [5] both with (green) and without (orange) channeling for DAMA. Channeling 
greatly reduces the scattering cross-sections necessary to produce the DAMA signal in the 
lower WIMP mass region (~9-15 GeV) and opens up even lower mass regions (down to 
~2 GeV), albeit only at a very low 7<r C.L. CDMS and XENON10 exclude all the 3a C.L. 
best fit regions. The low exposure experiments — CoGenT, CRESST I, and TEXONO — do 
not quite constrain the DAMA best fit regions. Both our older analyses and those performed 
by some of the experimental groups indicate these experiments do constrain portions of the 
DAMA region, but there are several differences between those analyses and our analysis here 
that change the DAMA best fit region (aside from any potential conservative assumptions 
on the experimental data sets as described in Section HTB : (1) an increase in the number 
of DAMA data bins from 2 to 36; (2) inclusion of DAMA's finite energy resolution, which 
increases the number of scatters over threshold and reduces the necessary cross-section, 
particularly at lower WIMP masses; (3) inclusion of the channeling effect, where indicated; 
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and (4) no background subtraction. Aside from shifting the DAMA best fit region downward, 
these changes lead to the very low mass regions disappearing from the best fit confidence 
regions, where the low threshold experiments provide the strongest constraints. SuperK 
provides relatively weak constraints for this SI case. 

The best fit SI cross-sections at low WIMP masses are indicated by the raster scan regions 
in Figure M (without channeling) and Figure [9] (with channeling). Without channeling, 
CRESST I and, to some degree, TEXONO constrain the 3a best fit cross-sections at masses 
below about 3 GeV. The best fit cross-section at these masses presents a poor fit to the 
DAMA data: they are incompatible at at least the 5a level. With channeling included, the 
best fit cross-sections shift downwards at these low masses and are no longer constrained by 
TEXONO. CRESST I excludes the best fit cross-sections only below ~1.5 GeV in this case. 
At high WIMP masses (above ~100 GeV), with and without channeling, the best fit cross- 
section is actually zero. The reason for this is that, for heavy WIMPs, the lowest energy 
bins will actually have a negative modulation amplitude, whereas the observed amplitude 
is positive. The best fit cross-section would then be an unphysical negative value; the best 
physically allowed value is zero. As at very low masses, the best fit cross-sections at these 
high masses represent poor fits to the DAMA data (incompatible at > 5a). 

To examine the compatibility of DAMA with other experiments, we use the x 2 g-o.f. 
regions shown in Figure [10] (without channeling) and Figure [TT] (with channeling). Without 
channeling, CDMS and XENON10 exclude all the parameter space that is compatible with 
the DAMA data within the 3a level. With channeling, CDMS and XENON10 exclude all 
the parameter space that is compatible with the DAMA data within the 90% level. There 
are parameters at 8-9 GeV that are compatible with DAMA within the 3a level and still 
satisfy the other constraints. Lower masses, while not constrained by the null experiments 
or DAMA's total rate, are only compatible with the DAMA data below the 3a level. 

We note that these results are sensitive to the thresholds of the two large exposure exper- 
iments, CDMS and XENON10. In particular, the DAMA 3a g.o.f. region with channeling 
would have been excluded with the use of the original £ e g value of 0.19 for XENON10, 
instead of the more recent measurement of 0.14 (which shifted the XENON10 threshold 
from 4.6 keV to 6.1 keV). The small window of compatibility at low WIMP masses is thus 
very dependent on the calibrations of these experiments' detectors. Given the closeness of 
the upward turn-off in the CDMS and XENON10 constraints (due to their thresholds) with 
the low mass DAMA region, even small adjustments to the calibrations, which will shift the 
turn-off in the constraints to slightly higher (or lower) WIMP masses, may open up more of 
the DAMA low mass region (or exclude it completely). 

F. Spin-Dependent Scattering 

We here present results for the case of purely spin-dependent interactions. We study the 
cases of scattering off of protons only, neutrons only, and mixed proton & neutron scattering. 

1. Proton- Only 

For the SD proton-only scattering case (a n = 0), the DAMA best fit parameters and null 
experiment exclusion limits are shown in Figure [61 Unlike the SI case (and the SD neutron- 
only case to be examined below), there is no high mass best fit region at better than a 5a 
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FIG. 12: Experimental constraints and DAMA best fit parameters for spin-dependent (SD) proton- 
only scattering. The DAMA raster scan and likelihood ratio regions are as described in Figure 
Null experiment constraints are as given in Figured! 
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FIG. 13: Same as Figure fT2| but including the channeling effect for DAMA. 
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FIG. 14: Experimental constraints and DAMA best fit parameters for spin-dependent (SD) proton- 
only scattering. The DAMA raster scan and x 2 g-o.f. regions are as described in Figure [TU1 Null 
experiment constraints are as given in Figure [6l 




FIG. 15: Same as Figure HU but including the channeling effect for DAMA. 
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C.L.; the best fits are strictly found around 11 GeV. The addition of channeling does not 
reduce the cross-sections of the best fit region as much as in other cases. These differences 
are a result of the relatively small contribution from scatters off of I in this case: where I 
presents the best fit to the DAMA data (the high mass region), there are enough Na scatters 
(which do not match the observed spectrum at those masses) to make the total spectrum 
a relatively poor fit, thus eliminating the higher mass region found in other cases at the 
90% & 3a C.L.'s. The relatively small amount of I scatters means less of the channeled I 
scatters, which contribute to the low mass region, so the low mass region does not shift to 
much lower cross-sections. 

Constraints from experiments with proton-even detector materials (e.g. Si, Ge, and Xe) 
are significantly weakened as those elements have very little spin in their proton groups for a 
WIMP to couple to. Noteably, CoGeNT and TEXONO (both Ge) constraints are too weak 
to appear in the various SD proton-only figures. CDMS and XENON only mildly constrain 
the 5<t and 7a C.L. DAMA best fit regions at higher WIMP masses. CRESST I, which uses 
proton-odd Al, still has reasonably strong constraints, but these fall just outside the 7a C.L. 
DAMA best fit region. SuperK, on the other hand, strongly constrains all of the DAMA 
best fit regions above their analysis threshold of 10 GeV. Only a small portion of the 3a 
C.L. region survives around 9-10 GeV, but would almost certainly be excluded if SuperK 
extended their analysis to lower masses. If any of the assumptions required for analyzing 
indirect detection limits from high energy solar neutrinos does not hold (see Section [Till) , 
so that the SuperK constraint does not apply, nearly all the the DAMA best fit parameter 
space survives in this case. 

Not shown in the figures are other proton-odd experiments such as COUPP [l(J 11] and 



KIMS [14| that can provide additional constraints in this SD proton-only scattering case. 



As noted in Section IHH COUPP does not present their results in a way that allows us 
to perform our own analysis and, hence, has not been included here. Both COUPP and 
KIMS would exclude the higher WIMP mass regions compatible with DAMA. However, 
KIMS has a relatively high threshold and cannot constrain the low WIMP mass DAMA 
region. COUPP could potentially constrain some of this low mass region, although without 
a reanalysis it is unclear to what degree (if it all) COUPP actually constrains the DAMA 
regions at these low WIMP masses. 

The best fit SD proton-only cross-sections at each mass are indicated by the raster scans 
in Figure [T21 (without channeling) and Figure [TBI (with channeling). CRESST I strongly 
constrains some of the raster scan band at low masses (<3 GeV) without channeling and 
only mildly constrains that band at low masses with channeling. Again, the best fit cross- 
sections at those low masses, as well as very high masses, are a poor fit to the data. 

Parameters compatible with the DAMA data are indicated by the x 2 g-o-f- regions shown 
in Figure [13] (without channeling) and Figure [151 (with channeling). Regions compatible with 
each experiment at the 90% level are found at 8-10 GeV both with and without channeling 
1 . In addition, regions compatible with DAMA within the 3a level can be found that satisfy 
the DAMA constraint at 7-10 GeV, or as low as 6-10 GeV when channeling is included. 
Without the SuperK constraint, parameters compatible with DAMA within the 90% level 
can be found over 8-15 GeV. Within the 3a level, that range expands to 7-19 & 35-55 GeV 



7 A cautionary note: saying parameters are consistent with each experiment at e.g. the 90% level is not 
equivalent to saying the experiments are consistent within the 90% level. The latter is a statement of a 
combined analysis, which is not performed here. 
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FIG. 16: Experimental constraints and DAMA best fit parameters for spin-dependent (SD) 
neutron-only scattering. The DAMA raster scan and likelihood ratio regions are as described 
in Figure O Null experiment constraints are as given in Figure [71 

(no channeling) or 6-18 GeV (with channeling), where some of the compatible regions have 
been limited by DAMA's total event rate. However, the high mass (35-55 GeV) DAMA 
compatible region without channeling would be incompatible with the COUPP and KIMS 
experimental constraints. 

2. Neutron- Only 

In the SD neutron-only coupling case (a p = 0), there is again two best fit mass regions, 
as shown in Figure [7J The addition of the channeling effect shifts the best fit cross-sections 
in the low mass region to lower values. For this neutron-only coupling case, experiments 
with the neutron-odd isotopes {e.g. Si-29, Ge-73, Xe-129, and Xe-131) are expected to gain 
in sensitivity relative to DAMA, which has only neutron-even isotopes (Na-23 and 1-127). 
The CDMS and XENON10 constraints exclude the ha C.L. DAMA best fit region with and 
without scattering. CoGeNT also constrains some of the 3a C.L. DAMA best fit region. 
Only DAMA best fit regions below a 5a C.L. survive the other experimental constraints. 
SuperK provides no constraints in this case. 

Once again, the DAMA best fit SD neutron-only cross-sections at each mass are indicated 
by the raster scans regions, shown in Figure [H] (without channeling) and Figure (T7J (with 
channeling). The best fit cross-sections within a 3<r C.L. are excluded by CRESST I for 
masses below 3 GeV and by CDMS and XENON10 for masses above 7 GeV; TEXONO 
constrains some of the best fit cross-sections at WIMP masses in between. For masses of 
3-7 GeV, the best fit cross-sections are consistent with the DAMA data only within the 
3a-5cr level. 

Parameters compatible with the DAMA data are indicated by the y 2 g.o.f. regions shown 



37 




M WIMP (GeV) 



FIG. 17: Same as Figure [16] but including the channeling effect for DAMA. 

in Figure [18] (without channeling) and Figure [19] (with channeling) . Parameters compati- 
ble with the DAMA data within the 3a level are excluded by the CDMS and XENON10 
constraints both with and without channeling. 

3. Mixed Couplings 

While the previous two cases took one of the two SD couplings to be zero, there is little 
reason to suspect one of the couplings to be significantly smaller than the other. Now we 
will consider the case of arbitrary SD couplings, where both couplings may be non-zero. 
We show in Figure [20] the experimental constraints and DAMA best fit parameters (WIMP 
mass and SD WIMP-proton cross-section) in the case of equal SD couplings (a n = a p ) and 
in Figure [21] for a n = — a p . In both cases, the SD scattering cross-sections off of the proton 
and neutron are equal. The null experiment constraints are nearly the same in both cases, 
but the DAMA best fit regions have some differences. In the compared 
to the the inclusion of channeling does not lower the best fit cross-sections in 

the low mass region by as much and the the higher mass region falls under a weaker CL. 
As with the previous cases, we find DAMA best fit regions within a 3a CL. survive other 
constraints at WIMP masses of ~8-10 GeV. Those regions fall within the 90% x 2 g-o.f. 
DAMA compatibility region (not shown), so there are parameters in these two cases that 
are compatible with each of the experiments within the 90% level. 

In Figure [22] we show the DAMA best fit couplings in the a p -a n plane, along with 
null experiment constraints, at several different WIMP masses; SuperK constraints have not 
been included. The gray regions indicate the regions allowed by the various null experiments, 
while the green and orange regions indicate the parameters most preferred (at a 3a CL.) for 
producing the DAMA signal with and without channeling, respectively (there is no DAMA 
best fit region at a WIMP mass of 5 GeV). The parameters in the overlapping regions (shown 
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FIG. 18: Experimental constraints and DAMA best fit parameters for spin-dependent (SD) 
neutron-only scattering. The DAMA raster scan and x 2 g.o.f. regions are as described in Fig- 
ure [TUJ. Null experiment constraints are as given in Figure [71 
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FIG. 19: Same as Figure [T8| but including the channeling effect for DAMA. 
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FIG. 20: Experimental constraints and DAMA best fit parameters for SD only scattering in the 
case a n = a p . The DAMA best fit regions are determined using the likelihood ratio method with 
(green) and without (orange) the channeling effect. 
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FIG. 21: Experimental constraints and DAMA best fit parameters for SD only scattering in the 
case a n = — a p . The DAMA best fit regions are determined using the likelihood ratio method with 
(green) and without (orange) the channeling effect. 
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FIG. 22: Experimental constraints and DAMA best fit parameters for SD only scattering in the 
a p -a Q plane at several different WIMP masses. The 3a best fit DAMA parameters are shown in 
orange (no channeling) and green (with channeling). The regions allowed at 90% by a combination 
of the null experiments, excluding Super-Kamiokande, is indicated by the gray regions. Parameters 
in the overlapping regions (gray+orange or gray+green) are consistent with both the DAMA best 
fit regions and all null experiments. Super-Kamiokande constraints, if valid, would exclude all of 
the overlapping regions shown at 12 GeV (and all such regions with WIMP masses above 10 GeV). 

in lighter shades of green and orange) are DAMA best fit parameters that satisfy the other 
constraints. Such regions can be found at masses of ~8— 15 GeV. Above 10 GeV, these 
regions have |a p | ^> |a n |. At masses of 8-10 GeV, |a p | > |a n |. SuperK constraints exclude all 
the overlapping regions above 10 GeV. 

Figure [23] similarly shows constraints in the a p -a n plane, but shows the couplings com- 
patible with the DAMA data within the 3a level using a \ 2 g-o-f- test. Again, SuperK 
constraints are not included. Overlapping regions in this case are the parameters that are 
compatible with all experiments (within 3a for DAMA and 90% for the others). Such regions 
are found over ~7-19 & 35-55 GeV without channeling and ~6-18 GeV with channeling. 
Inclusion of the SuperK constraint excludes all the overlapping regions above 10 GeV. 

G. Mixed SI & SD Scattering 

The SI only scattering case provides a slightly better fit to the DAMA modulation data 
than SD only scattering, as can be seen in Table IIV] although the improvement is mild 
enough that each of the SI and SD spectra provide reasonable fits to DAMA. In the case of 
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FIG. 23: Null experiment and DAMA constraints for SD only scattering in the a p -a n plane at 
several different WIMP masses. The 3a allowed DAMA parameters are shown in orange (no 
channeling) and green (with channeling). The regions allowed at 90% by a combination of the 
null experiments, excluding Super-Kamiokande, is indicated by the gray regions. Parameters in 
the overlapping regions (gray+orange or gray+green) are consistent with all experiments. There 
is no DAMA region at 5 GeV without channeling that is compatible within the 3a level. Super- 
Kamiokande constraints, if valid, would exclude all of the overlapping regions shown at 12 GeV 
and 60 GeV (and all such regions with WIMP masses above 10 GeV). 

no channeling, SI only scattering actually provides the best global fit, but mixed SI and SD 
couplings can also provide reasonable fits (only slightly worse than the SI only case). With 
channeling, the best fit has a mix of SI and SD contributions. 

Due to the linearity of the modulation amplitude and recoil events in each of the SI and 
SD cross-sections, a = ctsi + osd, one may find approximate mixed SI and SD couplings by 
interpolating between the SI and SD only regions given in the various figures. For example, 
at a WIMP mass of 10 GeV, an Si-only WIMP-proton cross-section of 1.5 x 10~ 5 pb and an 
SD proton-only WIMP-proton cross-section of 1 pb each fall within the 90% g.o.f. constraints 
of their respective cases. Then the case where ~60% of the signal comes from SI scattering, 
with cr Pi si = 0.9 x 10 -5 pb (60% of the cross-section that gives the full signal in SI), and 
the other ~40% comes from SD scattering, with ct Pi sd = 0.4 pb (40% of the cross-section 
that gives the full signal in SD), should fall within (or very near) a 90% g.o.f. constraint in 
the mixed SI & SD coupling case. A variety of mixed SI & SD couplings for WIMP masses 
around ~8 GeV can be found that reasonably fit the DAMA data as well as the various null 
results of other experiments. 
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VI. CONCLUSIONS 

The latest DAMA/LIBRA annual modulation data, which may be interpreted as a signal 
for the existence of weakly interacting dark matter (WIMPs) in our Galactic Halo, are 
examined in light of null results from other experiments. We investigate consequences for 
both spin-independent (SI) and spin-dependent (SD) WIMP couplings as well as mixed SI & 
SD couplings. The positive signal for annual modulation is compared to negative results from 
CDMS, CoGeNT, CRESST I, TEXONO, XENON10, and Super-Kamiokande (SuperK). In 
addition to the new experimental results, two new ingredients are included compared to 
our previous studies: 1) the energy spectrum of DAMA/LIBRA data is now given in 36 
bins, and 2) the effect of channeling is taken into account, which lowers the nuclear recoil 
energy threshold and increases sensitivity to low WIMP masses. Several statistical tools 
are implemented in our study: best fit regions are found; goodness-of-fit (g.o.f.); and raster 
scans over the WIMP mass. These approaches allow us to differentiate between the most 
preferred (best fit) and allowed (g.o.f.) parameters. It is hard to find WIMP masses and 
couplings consistent with all existing data sets; the surviving regions of parameter space are 
found here. For SI interactions, the best fit DAMA regions of any WIMP mass are ruled 
out to 3a, even with channeling taken into account. However, some parameters outside 
the best fit regions still yield a moderately reasonable fit to the DAMA data: we find 
parameters in the SI case at masses of ~8-9 GeV that are compatible with each of the 
experiments at better than the 3a level when channeling is included. For SD interactions, 
we find no neutron-only couplings that are compatible with the DAMA data to within the 3a 
level and satisfy other experimental constraints, with or without channeling included. For 
SD proton-only couplings, a range of masses below 10 GeV is compatible with DAMA and 
other experiments, with and without channeling, when SuperK indirect detection constraints 
are included; without the SuperK constraints, masses as high as ~20 GeV are compatible. 
Mixed couplings with a n = ±a p are found to be consistent with all experiments at ~8 GeV. 
The general case of mixed SD proton and neutron couplings are examined as well. In short, 
there are surviving regions at low mass for both SI and SD interactions (although only 
moderately so for the SI case) and at high masses for the case of SD proton-only couplings 
if indirect detection limits are relaxed. 

While in the process of preparing this paper, several papers were presented that also 
examine the compatibility of the DAMA signal with other negative experimental results, 



always assuming that the dark matter consists of WIMPs [26|, l67|, l68l. 169 



other non-WIMP candidates were studied as well, such as mirror [72j], composite [73[], and 



20, |21| (although 



WIMPless [74] dark matter). Our results are compatible with the results of these other 
studies of WIMPs, although our conclusions sometimes differ. In all of them channeling is 
taken into account. Ref. [26|] concludes that the recent CDMS results are compatible with the 
DAMA regions as derived by the DAMA collaboration, for WIMPs with spin independent 
interactions with mass in the 7 to 10 GeV range, a results similar to ours. Ref. [67J used 
only the raster scan and the DAMA data binned into two bins to study WIMPs with spin 
independent interactions (while here we use 36 data bins with this method) and find light 
WIMPs compatible with all experimental results, when channelling is included. Refs. [68[ 
and [70] both studied WIMPs with spin independent interactions using global fits with 
two parameters and the data binned into many bins, thus including the spectral modulation 
amplitude information. Their results are similar to those of our global fits, but they conclude 



that light WIMPs with spin independent interactions "cannot account for the data" [68 



or 



43 



"are strongly disfavoured" 1701 within the standard halo model (they considered also the 
limits imposed by the DAMA total spectrum, not included here, which, however, affect 
heavier WIMPs). Ref. [7l| studies indirect detection bounds which may be very important 
for light WIMPs, but depend on the WIMP-annihilation channels. 
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